/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */



/*! \file Wave.h
    \brief Multidimensional generalization of finite difference Laplacian operators.

	Multidimensional generalization of the classical Laplacian difference operator
	\f[
			\delta_{x\bar{x}} u_i=\frac{u_{i-1}+2u_i-u_{i-1}}{h^2}
	\f]
	The geometry is the unit cube \f$[0,1]^N\f$, where N is given by data::Dim
*/
#if !defined(_FD_Wave_)
#define _FD_Wave_
#include <src/miscellaneous/OS_basic.h>
#include <src/math/Operators/D_x.h>
#include <src/math/Operators/D__x.h>
#include <src/math/Operators/Delta.h>
#include <src/algorithm/Evolution_Algorithm/CFL_traits.h>
#include <src/algorithm/Evolution_Algorithm/Predictor_traits.h>

/*! 
    \brief Common multidimensional behavior for all finite difference linear p-Wave operators.

	Common multidimensional behavior for all finite difference Laplacian operators, i.e. for an input vector \f$\{v_i,w_i\}_i\f$ 
	return
	\f[
			\{\Nabla_x w_i,\Nabla_{\bar{x}} \cdot v_i\}_i
	\f]
	Where \f$F\f$ is a non linear operator. 
	The geometry is the unit cube \f$[0,1]^N\f$, where N is given by data::Dim
*/

template <class data, std::size_t = data::Dim>
class Wave 
	: public FiniteDifferenceOperator<data,data>
 {
//! return " \delta_{\bar{x}x}u"\f$.
		static std::string get_name(){return " /Wave_{/bar{x}x}u";};
		OS_STATIC_CHECK(false); // to implement
};



/*! 
    \brief One dimensional linear Wave operators Specialization.

	Common multidimensional behavior for all finite difference Laplacian operators, i.e. for an input vector \f$\{v_i,w_i\}_i\f$ 
	return
	\f[
			\{\delta_x w_i,\delta_{\bar{x}} v_i \}_i
	\f]
	The geometry is the unit cube \f$[0,1]\f$. 
*/
template <class data>
class Wave<data,2>  
	: public FiniteDifferenceOperator<data,data>
{
public :
		enum:std::size_t												{Order = 2};
		typedef OTSCentered												Centered_Type;
		typedef typename data::Boundary_type							Boundary_type;
		typedef typename data::value_type								value_type;
		
		typedef FiniteDifferenceOperator								Base;
		typedef typename Base::graph_type								graph_type;

		typedef  typename Base::Base_Cast								Base_Cast;
		typedef Wave													type;

		//! return " Wave u"\f$.
		static std::string get_name(){return " Wave u";};
		Wave( ) {};
		Wave(const Wave & data) : D_x_(from.D_x_), D__x_(from.D__x_), Base(data) {};
		Wave(const weight_set_type & data) :  Base(data) {};
		Wave(const Base_Cast &first) : Base(first)
		{};
		virtual ~Wave(){};
		Wave & operator = (const Wave & from) {
			D_x_ = from.D_x_;
			D__x_ = from.D__x_;
			Base::operator =(from);
			return *this;
		};
		virtual const graph_type& get_graph() const {
			return D__x_.get_graph();
		};
		virtual void update() {
			D_x_.set_state_variable(get_state_variable());
			D__x_.set_state_variable(get_state_variable());
		};
		virtual result_type operator()(typename argument_type from) const 
		{
			result_type result(from.size(2),from.size(1));
			result[0] = D_x_(from[1]);
			result[1] = D__x_(from[0]);
			return result;
		};
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
		};
		D_x<value_type, value_type::Dim, smart_ptr<data> >	D_x_;
		D__x<value_type, value_type::Dim, smart_ptr<data> >	D__x_;
};

template <class data>
class Wave_visc  
	: public Wave<data>
{
public :
		typedef Wave_visc									type;
		typedef Wave										Base;
		typedef typename Base::value_type					value_type;
		typedef typename Base::graph_type					graph_type;

		//! return " Wave u"\f$.
		static std::string get_name(){return " Wave_Viscos u";};
		Wave_visc( ) {};
		Wave_visc(const Wave_visc & data) : Base(from) {};
		Wave_visc(const weight_set_type & data) :  Base(data) {};
		Wave_visc(const Base_Cast &first) : Base(first)
		{};
		virtual ~Wave_visc(){};
		Wave_visc & operator = (const Wave & from) {
			Base::operator =(from);
			return *this;
		};
		virtual void update() {
			Base::update();
			Delta_.set_state_variable(get_state_variable());
		};
		virtual result_type operator()(typename argument_type from) const 
		{
			result_type result(from.size(2),from.size(1));
			result[0] = D_x_(from[1]);
			result[1] = D__x_(from[0])+Delta_(from[1])/(from.size(1));
			return result;
		};
		virtual const graph_type& get_graph() const {
			return Delta_.get_graph();
		};
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
		};
		Delta<value_type, value_type::Dim, smart_ptr<data> >	Delta_;
};

/*! \brief Bind operator of Dalembertian difference operator

	This is the time Bind operator of Laplacian difference operator Wave_visc, that is Wave_visc(t,u)=Wave_visc(u)
*/
template <class data,typename Time = OS_double >
struct Wave_visc_t : public BindTimeOperator<Wave_visc<data> ,Time > {

	typedef Wave_visc_t												time_operator_type;
	typedef BindTimeOperator<time_operator_type ,Time >				Base;
};


/*! \brief Bind operator of Laplacian difference operator Delta_traits

	This is the time Bind operator of Laplacian difference operator Delta_traits, that is Delta_t(t,u)=Delta(u)
*/
template <class data, typename F = Null, typename Time = OS_double >
struct Wave_t : public BindTimeOperator<Wave<data> ,Time > {

	typedef Wave_t														time_operator_type;
	typedef BindTimeOperator<time_operator_type ,Time >					Base;
	Wave_t(){};
};



#endif